Global behavior of cosmological dynamics with interacting Veneziano ghost 
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In this paper, we shall study the dynamical behavior of the universe accelerated by the so called 
Veneziano ghost dark energy component locally and globally by using the linearization and nullcline 
method developed in this paper. The energy density is generalized to be proportional to the Hawking 
temperature defined on the trapping horizon instead of Hubble horizon of the Friedmann-Robertson- 
Walker (FRW) universe. We also give a prediction of the fate of the universe and present the 
bifurcation phenomenon of the dynamical system of the universe. It seems that the universe could 
be dominated by dark energy at present in some region of the parameter space. 
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I. INTRODUCTION 

Currently, lots of theoretical models were proposed to 
explain the accelerating of the universe, which is con- 
vinced by many observations. In some of the models, 
people proposed a new kind of dark energy component 
with negative pressure in our universe, which will drive 
the acceleration, and the simplest dark energy model is 
the cosmological constant, but it suffers fine-tuning and 
coincidence problems. While, in other models, people try 
to modify the Einstein gravity at large scale in the uni- 
verse, e.g. f(R), DGP, etc. model, then the universe can 
be accelerated without introducing dark energy. 

Recently, a very interesting dark energy model called 
Veneziano ghost dark energy has been proposed [l[ , and 
in this model, one can obtain a cosmological constant 
of just the right magnitude to give the observed expan- 
sion from the contribution of the ghost fields, which are 
supposed to be present in the low-energy effective theory 
of QCD without introducing any new degrees of free- 
dom. The ghosts are needed to solve the U(l) prob- 
lem, but they are completely decoupled from the phys- 
ical sector 0. The ghosts make no contribution in the 
flat Minkowski space, but make a small energy density 
contribution to the vacuum energy due to the off-set of 
the cancellation of their contribution in curved space or 
time-dependent background. For example, in the Rindler 
space, the contribution of high frequency modes is sup- 
pressed by the factor e - 27Tk / a T anc [ the mam contribution 
comes from k ~ ax, where or is the temperature on the 
horizon seen by the Rindler observer 0] . In the cosmolog- 
ical context, one can choice ax ~ H , which corresponds 
to the temperature on the Hubble horizon. Then, in 
the context of strongly interacting confining QCD with 
topological nontrivial sector, this effect occurs only in 
the time direction and their wave function in other space 
directions is expected to have the size of QCD energy 
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scale. Thus, this ghost gives the vacuum energy density 
proportional to A^ CD H . With H ~ 10~ 33 eV, it gives 
the right order of observed magnitude ~ (3 x 10 _3 eV) 4 
of the energy density. 

It should be emphasized that the Veneziano ghost from 
the ghost dark energy model is not a new propagating 
physical degree of freedom and the description of dark 
energy in terms of the Veneziano ghost is just a matter 
of convenience to describe very complicated infrared dy- 
namics of strongly coupled QCD and it does not violate 
unitarity, causality, gauge invariance and other impor- 
tant features of renormalizable quantum field theory, see 

HE. 

Generally, it is very difficult to accept the linear be- 
havior that the energy of FRW universe is linear in Hub- 
ble constant "H", because QCD is a theory with a mass 
gap determined by the energy scale lOOMeV. So, it is 
generally expected that there should be an exponentially 
small corrections rather than that linear corrections H. 
However, as the arguments discussed in refs.[4|-|7| that 
the linear scaling H is due to the complicated topologi- 
cal structure of strongly coupled QCD, not related to the 
physical massive propagating degree of freedom. There- 
fore, the linear in Hubble constant H scaling has a strong 
theoretical support tested in a number of models. For 
recent paper on fitting the ghost dark energy model with 
current observational data including Snla, BAO, CMB, 
BBN, Hubble parameter and growth rate of matter per- 
turbation, see [8[ and also see Q. For recent progress on 
dark energy, see flol ]. 

Einstein equations can be written in a form called " uni- 
fied first law" based on the general definition of black 
hole dynamics on trapping horizon, which was proposed 
by Hayward There are some other horizons in the 
context of the FRW universe, such as the future inner 
trapping horizon defined in the next section, while the 
outer one is used to define black holes in a general space- 
time including time-dependent spacetime [111 ]. So, in this 
paper, we will choice ot ~ T t , where T t is the temper- 
ature defined on the trapping horizon of the FRW uni- 
verse instead to study its dynamical behavior locally and 
globally. In addition, we also study the bifurcation phe- 
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nomenon in this dynamical system and predict the fate 
of the universe. 

This paper is organized as follows: In Sec|lH we shall 
derive the equations of the autonomous system of the uni- 
verse with ghost dark energy component, and investigate 
its properties and solutions qualitatively and numerically 
in Sec lIIIl The global behavior of this system will be pre- 
sented in Sec lIVl and we also make a prediction of the 
fate of the universe in Sec|V] and study its bifurcation 
phenomenon in Sec lVIl In the last section, we will give 
some discussions and conclusions. 



II. AUTONOMOUS SYSTEM 

The spacetime of our universe is described by the FRW 
metric, which could be written in the form of 



ds 2 = h ab dx a dx b + rdfl 2 



(1) 



where x° — t, x 1 = r and r = a(t)r, which is the radius 
of the sphere while a(t) is the scale factor. Defining 



d£ = 7=\dtT ; 

V2 V VI - kr 2 



dr 



(2) 



where k is the spacial curvature, the metric could be 
rewritten as a double-null form 



ds 2 = -2d£ + dC + f- 2 dfl 2 , 



(3) 



then we get the trapping horizon np by solving the equa- 
tion d+f\f = f T = as 



, k 



-1/2 



(4) 



which coincides with the apparent horizon [12j . The sur- 
face gravity is given by 



K = -— { H + 2H 2 



1 



= --(l-e), (5) 



TT 



where 



2Hf T 



(lnfr)' 



(6) 



where dot and prime denote the derivative with respect 
to t and In a, respectively. Here, we assume e < 1 so that 
k < corresponding to inner trapping horizon. 

As we mentioned before, one can simply choice the 
temperature on the horizon as T ~ H in the cosmological 
context, but this is a special choice and generally, in the 
context of dynamical spacetime, the temperature is well- 
defined on the trapping horizon, which is proportional to 
the surface gravity on that horizon, namely T t — |k|/27t. 
Actually, it will reduce to T ~ H in the de Sitter space- 
time when we neglect the spatial curvature, because the 
trapping horizon is coincident with the Hubble horizon 



in this situation. In the following, we will also take the 
idea of the ghost dark energy model, in which the vac- 
uum energy density is proportional to the temperature 
defined on the trapping horizon pjje ~ T t . Therefore, 
the energy density of the Veneziano ghost dark energy is 
then given by 



Pde = ^ (1 - e) , 



(7) 



where C is a constant and here we have used Eq. ([5]). 
The Friedmann equation now reads 



k 



1 



H + —z = - [Pde + Pn 



(8) 



where p m is the energy density of dark matter, and by 
solving the Friedmann equation, we get 



Pde 



C 2 {l-e) 
6 




12 a, 



C 2 (l-e) 2 



(9) 



Usually, people also consider case that dark energy and 
dark matter are coupled to each other by ad ding a in- 
teraction term in the conservation equations [13] . Of 
course, the entirely independent evolution of the dark 
energy is a special case with a vanished interaction term. 
Furthermore, the microphysics seems to allow enough 
room for the coupling between the two components. At 
leat, before drawing a definitely conclusion, the inter- 
action case is still deserved to study. However, once 
a generic interaction term is introduced, the dynamical 
equations would become much more complicated and can 
not be solved analytically. As we known, by using the lo- 
cally dynamical analysis, one can get the future behavior 
of the nonlinear system and know whether the system is 
stable or not near the critical points. But, this method 
can only give the neighborhood properties of the criti- 
cal points. So, in the following, we will develop a global 
method to study the global behavior of the nonlinear dy- 
namical system, as well as its bifurcation phenomenon. 
As we known, the conservation equations including a gen- 
eral interaction between the dark energy and dark matter 
is described by 



Pde 



Pm + 3H Pm = Q, (10) 
3if(l + w DE )p DE = -Q, (11) 



where Q is the interacting term, and we will take its form 
as Q — 3H(apDE + Pp m ) in the following. Here a and (3 
are some constants. Usually, the existing literature only 
consider a simple case a = or /? = 0, as well as the 
case a = ft, see Refs. [l4[ and references therein. We 
will see that, once we have developed a global analysis 
method, it would be unnecessarily to set these special 
values, and there would be a interesting bifurcation phe- 
nomenon when one leaves these parameters free. 

The equation of state parameter of dark energy then 
be derived as 



1 

WDE = ~ . . 

3 \ 1 — e 



2e- 



Q 



PdeH 



- 1 . 



(12) 



3 



By using Eqs. 



TO]) and (fTTj) . we can obtain 

Pm) > (13) 



k If 

H ~ = --( pde(1 + wde) 

a z 2 \ 

and we also have the relation 

■ k _ 2e 

U J — ~— ■ 



(14) 



a r, 

By using the following set of dimensionless variables 

-k 



flDE — 



Pde 



n k = 



(15) 



then, the Friedmann equation ([8]) can be written as 

0-de + Mm + tt k = 1 , (16) 

which indicates that all terms in the above equation take 
values in the interval [0, 1]. The equation of state param- 
eter (fT2"]) could be rewritten in terms of new variables 
as 



wde 



3(1 



2e 



a - P[i - 1 . 



(17) 



where [i is the ratio of energy density of dark energy to 
dark matter, i.e. /i = flm/flnE, and the deceleration 
parameter q = —d/(H 2 a) is given by 

q = -l+n k +2e{l+n)£l DE = -(l-2e)(l+fi)n DE , (18) 



where we have used the constraint (|16|) to get the last 
equality of the above equation. By using Eq. (|13j) and 
(H3J), one can get the relation 



4e 



1 



W DE 
1 + M 



(19) 



Thus, the dynamical equations can be written in the fol- 
lowing form 



3 



4e\ 



(p + 1) 



4c 



2e + 3a + 3(/3-l + y)/x 



(20) 
(21) 



The Eqs. (120)) and (|2Tj) constitute an autonomous system 
and we shall study the properties of its solutions qualita- 
tively and numerically, for some recent works on locally 
dynamic analysis see Ref. [lj|, and also see [l6| for a 
review. 

Before ending this section, we would like to point out 
that the computation of the speed of sound should not 
be treated as a signal for instability of the theory as the 
Veneziano ghost is not a propagating physical degree of 
freedom. Such an "apparent signal for instability" is in 
fact a result of treatment of the Veneziano ghost as the 
conventional physical degree of freedom satisfying clas- 
sical equation of motion. Appropriate interpretation is 
that a notion of the speed of sound does not exist for a 
non-physical, non- propagating degree of freedom. This 
viewpoint was adapted in recent paper, see 



III. THE PICTURE OF TWO-DIMENSIONAL 
PARAMETER PLANE 

The critical points of the autonomous system can be 
obtained by setting // = and e' = 0. All the critical 
points are listed in coordinate (/j, c , e c ) in TabU And one 



CP. 


He 


e c 


P.C. (h > 0) 


Pi 


at 
1-/3 





a>0,/3<lora<0, /3>1 


p 2 


a 


1 


a < 0, P > -1/3 or a > 0,/3 < -1/3 


p 3 


-l 


3(q-/3 + 1) 


always unphysical 


Pa 


-l 


1 


always unphysical 



TABLE I: All the critical points (CP.) and their physical 
conditions (P.C). 

can see that if a > 0,/3 < 1 or a < 0,/3 > 1, then 
Pi is accepted from physical condition, namely, fi > 
by definition, and if a < 0,(3 > —1/3 or a > 0,/3 < 
— 1/3, then P% is accepted physically. To investigate the 
property of these critical points and the bifurcation of 
two-parameter family a and /3, we can write the variables 
near these critical points (/i c , e c ) in the form /i = /i c + S/i 
and e = e c +Se with the perturbations Sp and <5e, and thus 
we get the perturbation equations at the critical points 
Pi as follows 




= At 




(22) 



where Ai are the Jacobian matrices evaluated at the cor- 
responding critical point Pi, see Appendix [Ai for explicit 
forms of these matrices. 

The character equation of Ai can be written as 



A 2 - (trAj)A + det A L = 0, 



(23) 



and the location of (trA;, det Ai) relative to the parabola 
(tr Ai) 2 — 4 det Ai = in the trace-determinant plane de- 
termines the property of critical points Pi. A 2 x 2 matrix 
can be regarded as a 4-dimcnsional space determined by 4 
elements of the matrix. The trace-determinant plane is a 
2-dimcnsional representation of the 4-dimensional space. 
Thus, there are infinite different matrices corresponds to 
each point in the trace-determinant plane. Furthermore, 
each point in the a-f3 plane corresponds to some point 
in the trace-determinant plane, since Ai depends on the 
parameters a and j3. In other words, the property of 
critical points Pi may change with the parameters a and 
/3, because the family of dynamical system (|20| and (l2Tj) 
depends on a and /3. By using a lengthy and straight 
computing, we can divide the a-/3 plane into 11 domains 
(see FigfT]) where all the critical points have same type 
in each domain. In Tab. [TTJ all the critical points and 
their properties are listed. Noteworthily, the property of 
critical point is obtained from the linearized system of 
the nonlinear system near the critical point in the Tab. 
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III1 We say that a critical point Pj of a nonlinear sys- 
tem is hyperbolic if the eigenvalues of linearized system 
have nonzero real parts. The linearization theorem show 
that the nonlinear flow is conjugate to the flow of the 
linearized system in a neighborhood of Pj. Therefore, 
there properties are preserved for the nonlinear system 
and (p?Tj) . In the domains I and XI, there are two 
physical critical points P± and P%, and there is a physical 
point P 2 (Pi) in the domain III- VI (VI- IX). II and X are 
domains without physical points. 




-6 -4-2 2 



FIG. 1: The two-parameter plane is divided into 11 domains 
in accordance with the type of critical point and the physical 
condition /i < 0. All the critical points have same property 
in each domain. 



Domain 


Pi 


P 2 




Pa 


I 


source(p.c.p) 


source(p.c.p) 


saddle 


saddle 


II 


no physical critical point 


III 


sink 


saddle(p.c.p) 


saddle 


saddle 


IV 


spiral sink 


saddle(p.c.p) 


saddle 


saddle 


V 


saddle 


saddle(p.c.p) 


sink 


saddle 


VI 


saddle 


saddle(p.c.p) 


saddle 


sink 


VII 


saddle(p.c.p) 


saddle 


sink 


saddle 


VIII 


saddle(p.c.p) 


saddle 


saddle 


sink 


IX 


saddle(p.c.p) 


sink 


saddle 


saddle 


X 


no physical critical point 


XI 


saddle(p.c.p) 


sink(p.c.p) 


saddle 


saddle 



TABLE II: The critical points of the autonomous system (120[) 
and (|21[) and their type in each domain of the a- ft plane. 

In general, a bifurcation occurs when there is a re- 
markably change in the structure of the solutions of the 
dynamical system as the parameters a and ft. The sim- 
plest types of bifurcation occurs when the property of 
equilibrium solutions changes as a and ft vary. Thus, a 



bifurcation may be appearance when the parameters vary 
from some domain into neighbor one. For example, the 
structure of P1-P2 changes from saddle-saddle to saddle- 
sink and the structure of P3-P4 changes from sink-saddle 
to saddle-saddle when the location of (a, (3) crosses the 
line 3cv — 3/3 — 1 = 0, namely, (a, ft) moves from domain 
VIII to IX in the a- ft plane. 



IV. GLOBAL BEHAVIOR OF THE NONLINEAR 
SYSTEM 

In the section IIII1 we have used the technique of lin- 
earization to determine the behavior of solutions near 
physical points. In this section, we use the qualitative 
technique of nullcline for analyzing the global behavior 
of nonlinear system, which is one of the most useful tools 
for analyzing solutions of nonlinear planar system. For 
our system, the /i-nullclines are the set of points deter- 
mined by 

M = -i; (p-i + y) ? = -<*■ ( 24 ) 

Similarly, the e-nullclines are the set of points determined 
by 

e = l; + =~08-l)-«- (25) 

It is worth noting that the nullclines consist of lines and 
hyperbolas. The asymptotes of hyperbola in Eq. (pM|) are 
fi = and e = 1 — ft, which is interesting to determine 
the fate of universe. 

On the /i-nullclines, we have \J = 0, the vector field 
points (//, e') are vertical, so \i- nullclines divide /x-e plane 
into regions where the vector field points either to the left 
or to the right. Similarly, the e-nullclines separate the 
plane into regions where the vector field points either 
upward or downward. If we determine all of the null- 
clines, then this allows us to decompose fi-e plane into a 
collection of basic regions. The vector field must point 
in one of 4 directions (northeast, northwest, southeast or 
southwest) for any of the basic regions between the null- 
clines, since it is neither vertical nor horizontal. In each 
of basic regions, the vector field points are provided with 
a certain direction. Therefore, the basic regions allow 
us to understand the phase portrait from a qualitative 
point of view. Furthermore, the nullclines are depending 
on the parameters a and ft in our dynamical system, so 
we have to discuss each of domains I- XI separately. 

Here, we discuss domain I and XI in detail, since there 
are two physical critical points Pi and Pi in these do- 
mains. In the case of a, ft G Domain I, the nullclines 
divide fi-e plane into 13 basic regions in FigO By taking 
one point in each of these regions firstly and then decid- 
ing the direction of the vector field at that point, we can 
determine the direction of the vector field at all points in 
the basic region. Similarly, the nullclines also divide /i-e 



5 



plane into 13 basic regions for the case of a, j3 G Domain 
XI, see Fig|31 and we can decide the direction of the vec- 
tor field at all points in the basic region. Especially, we 
can find out the approximate behavior of solutions every- 
where in the plane, which has interest for the discussion 
of physical significance. For example, we will know that 
the matter component is a dominant position at late evo- 
lution if any solution tends toward infinity in the basic 
region marked A(B) in Figj3] with the northeast (south- 
east) direction. In other words, we will know the fate of 
universe and its initial conditions. 



By using above analysis and the numerical calculation, 
we can show the global behavior of the nonlinear system 
exactly. In FigJH we plot evolution trajectories in the /i-e 
plane where we have chosen the parameters a, P € Do- 
main XI. We also plot evolution trajectories when a, (3 G 
Domain I in FigjSJ It is worth noting that the physical 
admissible range is the only right-half plane (/t > 0), so 
we must choose initial conditions at the right-half plane. 
The detailed discussion of the physical arguments will 
appear in Section IVl 




FIG. 2: The fi-e plane are divided into 13 basic regions by 
the nullclines as a, f3 £ Domain I. The direction of the vector 
field at all points in the basic region are determined. 




FIG. 3: The fi-e plane are divided into 13 basic regions by the 
nullclines as a, f3 g Domain XL We can decide the direction 
of the vector field at all points in the basic region. 



V. FATE OF UNIVERSE 

From FigHl one can see that if a, /3 belongs domain 
I, and the initial values of /i and e satisfy the following 
condition 

P-l + lj^Ho + a > 0, (26) 
3^-1 + ^/^ + 260 + 30 > 0, (27) 

i.e. the initial point (/xq, eo) belongs to A (J B (A denotes 
the yellow region, while B denotes the green region), the 
universe will evolve to matter dominate finally, which 
suggests that the expansion will asymptotically come to 
a halt. If a,j3 belongs domain I, and the initial value of 
fi and e satisfy the following condition 



P - 1 + 



4e r 



4e, 



/to + a < . 



(28) 



3 1/3- l + ^J/ t0 + 2e + 3a < 0, (29) 

the dynamical evolution will be in unphysical region, so 
such kind of initial values are excluded by the physical 
condition fi > 0. 

From Fig|5j one can see that if a, j3 belongs domain XI, 
and if not only the initial value of satisfies the physical 
condition /i > 0, but also satisfies eo > 1, then, the sys- 
tem will evolute to Pi finally, namely, the fate of universe 
will be 



Si r, 



n 



DE 



(a>0,P<-h. (30) 



If we assumed that we have get the above limit at present, 
then, we can choose a and j3 as 



DE,0 



m,0 



(31) 



where D,de,o and ri m ,o denote the present values of en- 
ergy density of dark energy and matter respectively. If 
a, P belongs to domain XI, and the initial values satisfy 
e < with physical condition /i > 0, the dynamical sys- 
tem will tend to be /t — > 0, namely, the universe will 
be dominated by dark energy, and will be accelerating 
eternally. 
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FIG. 4: Choosing a, ft in domain XI, the evolution trajectories 
in the fi-e plane are plotted. 
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P 2 (-3a/4, 1), P 3 (-l,3a/2) and P 4 (-l,l) are all saddle 
points, and the physical discussion is not change except 
losing physical point Pi . In FigjSJ we plot evolution tra- 
jectories in the fi-e plane for a — —1.485 and (3 — 1. For 
the parameters a, /3 G Domain III, we also plot evolution 
trajectories in the fj,-e plane for a = —4 and j3 = 0.217 in 
Fig J7] Obviously, Pi changes into a sink from a source as 
the parameters a, f3 run in Domain II from Domain I, and 
Pi also becomes a non-physical point. Therefore, one can 
say that this is a typical bifurcation for the nonlinear sys- 
tem. However, the physical arguments of evolution have 
not remarkable change. 



_ -—^^v 




— ^ 
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FIG. 6: As a = —1.485 and = 1, we plot the evolution 
trajectories in the fi-e plane. 




FIG. 5: We plot the evolution trajectories in the fi-e plane as 
a, /3 are in domain XI. 



VI. BIFURCATION 




The bifurcation phenomenon occurs when there is a 
remarkable change in the structure of the solutions of 
the dynamical system as the parameter a and /3. Here, 
we consider that the structure of the solutions makes a 
change as the parameter j3 crosses line (9 = 1. 

Now, we consider the case a < and (8 = 1. The 
system is simplified to 




p! = (4eju + 3a)0 + l), 

e = (2e + 3a + 4e^)(l-e). 



(32) 
(33) 



FIG. 7: Fixing a — — 4 and /3 = 0.217, we plot the evolution 
trajectories in the fi-e plane. 



In this case, there are only three critical points: 
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VII. CONCLUSION 

In conclusion, we have studied the dynamical behavior 
of the universe accelerated by the generalized Veneziano 
ghost dark energy component locally and globally. We 
have found that in this system, there are four critical 
point but only two of them have physical meaning when 
the parameters are chosen in some proper regions to sat- 
isfy the physical conditions. We have used the technique 
of linearization to determine the behavior of solutions 
near physical points and used the qualitative technique 
of nullcline for analyzing the global behavior of this non- 
linear system, which is one of the most useful tools for 
analyzing solutions of nonlinear planar system. We have 
shown that the universe could be dominated by the dark 
energy at present if we choice a set of suitable parame- 
ters. We also give an example to show the bifurcation 
phenomenon in this interesting dynamical system. Be- 
fore ending the paper, we would like to emphasis that the 
qualitative technique of nullcline we developed in this pa- 
per is very powerful and could be used in any nonlinear 
dynamical system, especially in the planar system, so it 
deserves further studying. 
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Appendix A: Jacobian matrixes 



f-3(l + a-P) 



An 



\ 3(13-1) 



-, o . O o 12ct(3oi-3/3-l) 

1 - 3a + 3/3 l 1+30) p 2 



A, = 



(Al) 



(A2) 



-3(1 + a- f3) 



-§(l + 3a-3/3)(l + 2a-/3) 1 + 3a - 3^ 



(A3) 
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